%% Reproduce_Simulations.m
%  with this file it is possible to generate the simulations with each of
%  the parameters' set. It will generate a folder for each model. The
%  generation 1 models are the models of session1, the generation 2 models
%  are the models of session2_sham and the generation 3 models are the
%  models of session2_tbs.
%
% Estimated run time: 30 minutes with 4 cores.


%% INITIALIZATION and CHECK IF ALL NECESSARY FILES ARE PRESENT

clear
close all
clc

if (  exist('DataSet.mat','file')~=2 ...
   || exist('_egidio2iPurkinje.dat','file')~=2 ...
   || exist('_egidio2iGranule.dat','file')~=2 ...
   || exist('_egidio2iGolgi.dat','file')~=2 ...
   || exist('_egidio2iPurkinje.cfg','file')~=2 ...
   || exist('_egidio2iGranule.cfg','file')~=2 ...
   || exist('_egidio2iGolgi.cfg','file')~=2 ...
   || exist('EBCC_s1.exe','file')~=2 ...
   || exist('EBCC_s2.exe','file')~=2 )
    disp('ERROR: There are some missing files')
    return
end


%% CREATE THE WEIGHTS AND THE PARAMETERS FILE FOR SESSION1
load DataSet.mat
TrialDuration=0.8;
nrep=70;
CR_event=ones(3,44,nrep)*NaN;
CR_event_exp = ones(3,22,nrep)*NaN;

CR1=[];
CR2=[];
CR3=[];


Unici_sani=unique(Genio_Migliori_sani','rows');
for i=1:size(Unici_sani,1)
    UniciN_sani(i)=sum(ismember(Genio_Migliori_sani',Unici_sani(i,:),'rows'));
end
Numerosita_sani=UniciN_sani';
Unici_sani(:,10)=[];

% Generate parameters file
save file_param.dat Unici_sani  -ascii

% Generate initial weights
Unici_sani=unique(Genio_Migliori_sani','rows');
for i=1:size(Unici_sani,1)
    peso(1,1:2)=[8000 0.7];
    peso(2,1:2)=[38438 Unici_sani(i,7)];
    peso(3,1:2)=[24 0];
    peso(4,1:2)=[1200 Unici_sani(i,8)];
    peso(5,1:2)=[24 Unici_sani(i,9)];
    fid=sprintf('WeightsEBCC_%i.dat',i-1);
    filena=fopen(fid,'w');
    for k=1:5
        fprintf(filena,'%d %1.10f\n', peso(k,1), peso(k,2));
    end
    fclose(filena);  
end

% SIMULATE SESSION1 MODELS WITH WEIGHTS GENERATION EACH TRIAL
gen_num=1;
disp(['%%%% GENERATION NUMBER ' num2str(gen_num) ' %%%%'])
parfor chr_num=0:size(Unici_sani,1)-1
    %%% RUN EDLUT FOR GENERATION gen_num
    command = ['EBCC_s1',' ',num2str(gen_num),' ',num2str(chr_num),' 1'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end   
end

for chromo=1:size(Unici_sani,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    load CR.dat
    for i=1:nrep
        CR_ordered=CR(find(CR<=(i*TrialDuration)));
        CRS=CR_ordered(find(CR_ordered>(i*TrialDuration-TrialDuration)));
        CR_event(gen_num,chromo,i)=size(CRS,1);
        clear CRS CR_ordered
    end
    a(1,:)=CR_event(gen_num,chromo,:);
    CR1=[CR1;[repmat(a',1,Numerosita_sani(chromo))]'];
    clear CR i a
    cd('..')
end
CR_event_exp(gen_num,:,:)=[Spre,Tpre]';

% SIMULATE
parfor chr_num=0:size(Unici_sani,1)-1
    %%% RUN EDLUT FOR GENERATION gen_num
    command = ['EBCC_s1',' ',num2str(gen_num),' ',num2str(chr_num),' 0'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end   
end



%%  CREATE THE PARAMETERS FILE FOR SESSION2_SHAM
clear ans fid filena gen_num i k peso
Unici_sani2=unique(Genio_Migliori_sani2','rows');
for i=1:size(Unici_sani2,1)
    UniciN_sani2(i)=sum(ismember(Genio_Migliori_sani2',Unici_sani2(i,:),'rows'));
end
Numerosita_sani2=UniciN_sani2';
Unici_sani2(:,7)=[];

save file_param.dat Unici_sani2  -ascii

% SIMULATE SESSION2_SHAM MODELS WITH WEIGHTS GENERATION EACH TRIAL
gen_num=2;
disp(['%%%% GENERATION NUMBER ' num2str(gen_num) ' %%%%'])
parfor chr_num=0:size(Unici_sani2,1)-1
    %%% RUN EDLUT FOR GENERATION gen_num
    command = ['EBCC_s2',' ',num2str(gen_num),' ',num2str(chr_num),' 1'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end
end

for chromo=1:size(Unici_sani2,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    load CR.dat
    for i=1:nrep
        CR_ordered=CR(find(CR<=(i*TrialDuration)));
        CRS=CR_ordered(find(CR_ordered>(i*TrialDuration-TrialDuration)));
        CR_event(gen_num,chromo,i)=size(CRS,1);
        clear CRS CR_ordered
    end
    a(1,:)=CR_event(gen_num,chromo,:);
    CR2=[CR2;[repmat(a',1,Numerosita_sani2(chromo))]'];
    clear CR i a
    cd('..')
end
CR_event_exp(gen_num,1:11,:)=[Spost]';

% SIMULATE
parfor chr_num=0:size(Unici_sani2,1)-1
    command = ['EBCC_s2',' ',num2str(gen_num),' ',num2str(chr_num),' 0'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end
end

%%  CREATE THE PARAMETERS FILE FOR SESSION2_TMS
clear gen_num i
Unici_tms2=unique(Genio_Migliori_tms2','rows');
for i=1:size(Unici_tms2,1)
    UniciN_tms2(i)=sum(ismember(Genio_Migliori_tms2',Unici_tms2(i,:),'rows'));
end
Numerosita_tms2=UniciN_tms2';
Unici_tms2(:,7)=[];

save file_param.dat Unici_tms2  -ascii

% SIMULATE SESSION2_TMS MODELS WITH WEIGHTS GENERATION EACH TRIAL
gen_num=3;
disp(['%%%% GENERATION NUMBER ' num2str(gen_num) ' %%%%'])
parfor chr_num=0:size(Unici_tms2,1)-1
    %%% RUN EDLUT FOR GENERATION gen_num
    command = ['EBCC_s2',' ',num2str(gen_num),' ',num2str(chr_num),' 1'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end
end

for chromo=1:size(Unici_tms2,1)
    cd(['.\SAVED_FILES_Generation_' num2str(gen_num) '_Chromosome_' num2str(chromo-1)])
    load CR.dat
    for i=1:nrep
        CR_ordered=CR(find(CR<=(i*TrialDuration)));
        CRS=CR_ordered(find(CR_ordered>(i*TrialDuration-TrialDuration)));
        CR_event(gen_num,chromo,i)=size(CRS,1);
        clear CRS CR_ordered
    end   
    a(1,:)=CR_event(gen_num,chromo,:);
    CR3=[CR3;[repmat(a',1,Numerosita_tms2(chromo))]'];
    clear CR i a
    cd('..')
end
CR_event_exp(gen_num,1:11,:)=[Tpost]';

% SIMULATE
parfor chr_num=0:size(Unici_tms2,1)-1
    command = ['EBCC_s2',' ',num2str(gen_num),' ',num2str(chr_num),' 0'];
    status = system([command]);
    while status ~= 0             % Status is 0 when execution is finished and is ok!
        % While the execution of EDLUT is running, wait
    end
end

%% REPRODUCE TABLE II - SENSITIVITY AND SPECIFICITY VALUES

% MODELS VS EXPERIMENTAL MEDIAN
CR_event=ones(3,753,nrep)*NaN;
CR_event(1,1:417,:)=CR1;
CR_event(2,1:588,:)=CR2;
CR_event(3,1:753,:)=CR3;

TP=zeros(3,1);
TN=zeros(3,1);
FP=zeros(3,1);
FN=zeros(3,1);

gen_num = 1;
for chromo=1:417
    for i=1:nrep
        if CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,:,i))
            TP(gen_num)=TP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,:,i))
            TN(gen_num)=TN(gen_num)+1;
        elseif CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,:,i))
            FP(gen_num)=FP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,:,i))
            FN(gen_num)=FN(gen_num)+1;
        end
    end
end

gen_num = 2;
for chromo=1:588
    for i=1:nrep
        if CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,1:11,i))
            TP(gen_num)=TP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,1:11,i))
            TN(gen_num)=TN(gen_num)+1;
        elseif CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,1:11,i))
            FP(gen_num)=FP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,1:11,i))
            FN(gen_num)=FN(gen_num)+1;
        end
    end
end

gen_num = 3;
for chromo=1:753
    for i=1:nrep
        if CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,1:11,i))
            TP(gen_num)=TP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,1:11,i))
            TN(gen_num)=TN(gen_num)+1;
        elseif CR_event(gen_num,chromo,i) & ~median(CR_event_exp(gen_num,1:11,i))
            FP(gen_num)=FP(gen_num)+1;
        elseif ~CR_event(gen_num,chromo,i) & median(CR_event_exp(gen_num,1:11,i))
            FN(gen_num)=FN(gen_num)+1;   
        end
    end
end

% PERCENTAGES:

SENSITIVITY_1=(100*TP./(TP+FN))';
SPECIFICITY_1=(100*TN./(TN+FP))';

% EXPERIMENTAL SUBJECTS VS EXPERIMENTAL MEDIAN

TP=zeros(3,1);
TN=zeros(3,1);
FP=zeros(3,1);
FN=zeros(3,1);

gen_num = 1;
for chromo=1:22
    for i=1:nrep
        for exp=1:22
            if CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                TP(gen_num)=TP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                TN(gen_num)=TN(gen_num)+1;
            elseif CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                FP(gen_num)=FP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                FN(gen_num)=FN(gen_num)+1;
            end
        end
    end
end

gen_num = 2;
for chromo=1:11

    for i=1:nrep
        for exp=1:11
            if CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                TP(gen_num)=TP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                TN(gen_num)=TN(gen_num)+1;
            elseif CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                FP(gen_num)=FP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                FN(gen_num)=FN(gen_num)+1;
            end  
        end
    end

end

gen_num = 3;
for chromo=1:11
    for i=1:nrep
        for exp=1:11
            if CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                TP(gen_num)=TP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                TN(gen_num)=TN(gen_num)+1;
            elseif CR_event_exp(gen_num,chromo,i) & ~CR_event_exp(gen_num,exp,i)
                FP(gen_num)=FP(gen_num)+1;
            elseif ~CR_event_exp(gen_num,chromo,i) & CR_event_exp(gen_num,exp,i)
                FN(gen_num)=FN(gen_num)+1;
            end    
        end
    end
end

% PERCENTAGES

SENSITIVITY_2=(100*TP./(TP+FN))';

SPECIFICITY_2=(100*TN./(TN+FP))';
clc
disp('REPRODUCE TABLE II - SENSITIVITY AND SPECIFICITY VALUES')
disp('Model vs Experimental Median');
disp(['Sensitivity: ' num2str(round(SENSITIVITY_1))])
disp(['Specificity: ' num2str(round(SPECIFICITY_1))])
disp('Experimental Subjects vs Experimental Median');
disp(['Sensitivity: ' num2str(round(SENSITIVITY_2))])
disp(['Specificity: ' num2str(round(SPECIFICITY_2))])


